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Abstract 

Direct numerical simulations (DNS) and large-eddy simulations (LES) simulations of flow through a 
periodic channel with a constriction are performed using the dynamic Smagorinsky model at two 
Reynolds numbers of 2800 and 10595. The LES equations are solved using higher order compact 
schemes. DNS are performed for the lower Reynolds number case using a fine grid and the data are 
used to validate the LES results obtained with a coarse and a medium size grid. LES simulations are 
also performed for the higher Reynolds number case using a coarse and a medium size grid. The 
results are compared with an existing reference data set. The DNS and LES results agreed well with 
the reference data. Reynolds stresses, sub-grid eddy viscosity, and the budgets for the turbulent kinetic 
energy are also presented. It is found that the turbulent fluctuations in the normal and spanwise 
directions have the same magnitude. The turbulent kinetic energy budget shows that the production 
peaks near the separation point region and the production to dissipation ratio is very high on the order 
of five in this region. It is also observed that the production is balanced by the advection, diffusion, and 
dissipation in the shear layer region. The dominant term is the turbulent diffusion that is about two 
times the molecular dissipation. 


I. Introduction 

T he major objectives of any CFD code are to predict global aerodynamic quantities such as lift and drag as well 
as to predict local flow features such as skin friction, heat transfer, and pressure oscillations, accurately and 
efficiently for a wide range of problems. The flow physics are fundamentally governed by the unsteady three- 
dimensional Navier-Stokes (N-S) equations that state the three conservation laws for mass, momentum, and energy. 
Turbulent flow is characterized by the existence of a vast range of length and time scales ranging from the smallest, 
the Kolmogorov scale, to the largest, determined by the geometry. 1 2 The required computer resources and time 
constraints hinder the solution of these equations numerically for turbulent flows at high Reynolds numbers. 
Analysis 2 has shown that the number of grid points required increases as Re 9 4 and the number of time steps required 
increases as Re’ /4 . These challenging requirements restrict the solution of the full N-S equations to simple 
geometries such as channels and flat plates at low Reynolds numbers. The first successful direct numerical 
simulation 3 (DNS) was performed for a channel flow at a Reynolds number of 3000. With increasing computational 
capability, the application of DNS to more complex problems at high Reynolds numbers is currently being pursued 
by many researchers. The recent status of DNS in turbulent flow is reviewed by Moin and Mahesh. 4 However, DNS 
usage is limited to understanding the turbulent physics in canonical problems and to extending the findings for 
applications to turbulence model development. Hence, the current state of the art is to solve some approximate 
versions of the N-S equations that require modeling. 

The simplest and the most popular approximate method is the set of Reynolds- Averaged Navier-Stokes (RANS) 
equations, where the equations are derived for time-averaged quantities. The unclosed Reynolds stresses are 
modeled to close the equations. Wilcox 3 gives a good account of different models and their applications. In general, 
existing models provide good results compared to experiments in attached flows where turbulence is in equilibrium. 
However, prediction of separated flows with existing RANS turbulence models is still not very satisfactory. It is 
difficult to accurately capture both the separation and reattachment points (i.e., the size of the separation bubble). In 
particular, it is believed that the flow field dynamics near the separation point are not predicted correctly. The RANS 
computations performed for the flow over a two-dimensional hump, NASA Langley CFDVAL2004 Case 3 6 ’ 7 , 
consistently over-predicted the reattachment point compared to the experiment. The computed Reynolds stresses 
were 2 to 3 times smaller than the measured values. Computations performed 8 for two other massively-separated 
flow cases (the 2-D periodic hill 9 and the 3-D Ahmed body 10 ) using different turbulence models also predicted too- 
long separation bubbles. These computations further exemplify the deficiencies in predicting massively-separated 
flows using existing RANS models. 
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The second approximate method uses the Large-Eddy-Simulation (LES) equations, which are obtained for the 
large-scale motions by fdtering small scales. The unknown stresses introduced by the small scales (subgrid-scale 
stresses, or SGS) are modeled to close the equations. Sagaut 11 and Gamier et al. 12 give account of different SGS 
models and their applications in incompressible and compressible flows, respectively. This approach is 
computationally more expensive than the RANS approach, but yields more accurate results without many ad-hoc 
modifications to the SGS models. This methodology is being applied to a wide range of problems. Recently, several 
DNS and LES simulations have been performed for separated flows over two- and three-dimensional humps in 
channel flows. 13 " 20 

Flow through a periodic channel with a constriction (separated flow over a 2-D hill, Fig. 1) is widely used as a 
test case to investigate the flow physics of turbulent separated flows over smooth surfaces and to validate different 
RANS and LES solutions obtained on coarser grids with highly resolved DNS and LES data sets. Breuer et al. 19 
gave a brief history behind the selection of this geometry as a benchmark test case. They also have presented highly 
resolved LES simulation results for a range of Reynolds numbers 700 to 10595. Frohlich et al. 15 also performed a 
highly resolved LES at a Reynolds number of 10595. The results include mean flow profdes, Reynolds stress 
profiles, energy balances, and the instantaneous flow fields. Ziefle et al. ls performed LES simulations with a very 
coarse grid using the approximate-deconvolution SGS model for compressible flow at a Reynolds number of 2800. 
It was found that the grid resolution near the separation point is critical in capturing the reattachment point correctly 
in the simulations. Even though this geometry is attractive for numerical simulations because of the simple geometry 
and periodic boundary conditions in the streamwise direction, the flow becomes very unsteady and complex due to 
the varying inflow conditions at the crest that is continuously influenced by the previous hill. In the previous 
simulations, it took about 50-60 flow through times to obtain reasonable statistical averages. Another observation 
from the LES and DNS was that the ratio of production over dissipation takes very high values in the separated 
shear region compared to homogeneous flows. Whether this is one of the causes for the deficiency in RANS is not 
clear. The turbulent dynamics between the separation point and the reattachment point is very complex. 21 The 
adverse pressure gradient near the reattachment point and spreading of the shear layer due to turbulence are 
interconnected and the balance between these two factors determines the separation and the reattachmcnt points. 

The objective of the present research is to first validate our DNS/LES results for the periodic channel flow with a 
constriction with the reference data of Breuer et al. 17 ’ 19 Another objective is to investigate the turbulent dynamics of 
this flow by analyzing the Reynolds stresses and the budget of different terms in the turbulent kinetic energy 
equation. We also want to investigate the effects of the SGS model on the turbulent statistical quantities by 
performing the LES with and without the model on the same grid. We solve the three-dimensional unsteady 
compressible Navier-Stokes equations using higher-order compact schemes for space discretization and 3 rd order 
Runge-Kutta scheme for time integration. We employ the dynamic Smagorinsky SGS model in the LES 
calculations. We first present the DNS results obtained on a fine grid at a low Reynolds number based on the bulk 
velocity (£4), density (p 6 ), the hill height (h), and molecular viscosity at the wall (p,,,) of Ret = Pb Ut h/p w = 2800, 
and then present the wall-resolved LES (WRLES) results at low to high Reynolds numbers of Rei,= 2800 and 10595. 

I. Models and Flow Conditions 

Computations were performed for a flow through a channel with a constriction (2D-hill). The geometry and the 
coordinate systems are shown in Figs. 1(a) and 1(b). The flow was periodic in the streamwise and spanwise 
directions. The lengths are non-dimensionalized by the hill height, h. The computational dimension in the 
streamwise direction is L x = 9, in the spanwise direction is L : = 4.5, and in the normal direction above the crest of the 
hill is L y = 3.035. These are the same dimensions used in the DNS simulations 15 ' 19 and in the LES simulations. 15 
Simulations are performed for Reynolds numbers based on the bulk velocity, bulk density, and the hill height of Re* 
= pb Ub h/p w = 2800 and 10595. The wall temperature is fixed at a constant temperature of 300K. The Mach number 
based on the bulk velocity is M 0 = 0.2. The non-dimensional bulk density and the bulk velocity are defined as 
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Figure 1. Periodic channel with a constriction and the coordinate systems used. 


II. Governing Equations 

The partial differential equations solved are the filtered three-dimensional unsteady compressible Navier- 
Stokes equations in conservation form 22 
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Here the over-bar denotes filtered variables and the tilde denotes the Favre-averaged variables. The subgrid- 
scale (SGS) terms appearing in the momentum equations are 


3 


( 8 ) 
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Here xy is the subgrid-scale stress tensor. The subgrid-scale terms appearing in the energy equations are 
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The subgrid-scale viscous term Zy is orders of magnitude smaller than the subgrid-scale stress term and is neglected. 
The deviatoric part of the SGS stress tensor is modeled using the Smagorinsky eddy viscosity model 23 


T u -- T kAj = 


y y 

2 I 


v t =(QA ) 2 \S(u)\, 
|5(M)| 2 =i5,(«)5 /7 (M). 


(17) 


(18) 


The value of the parameter C s varies from 0.18 for isotropic turbulence to 0.10 for plane channel flow. 24 ' 23 The 
characteristic length scale A is determined by the local grid size. 


A = (AxAyAz) 


1/3 


(19) 


The model is very dissipative in transitional flows and does not vanish near the wall and in laminar shear flows. 
One has to use Germano’s 26 dynamic procedure or a simple Van Driest damping to remedy these issues. We 
followed Vreman’s 27 procedure in the derivation of the dynamic Smagorinsky model (DSM). In the dynamic eddy- 
viscosity model, the deviotoric stress is written as a linear function of a variable coefficient Cd. 
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(20) 


T°-C,A J |S(5)|S„(ii). 

The coefficient Cd is determined dynamically based on Germano’s identity, which relates the flow fields at the 
filtered level and at a new test-filter level. Applying the Germano identity to the subgrid-scale stress term yields the 
following system of equations for the coefficient Cd. 


where 


C d M ir L ip 

Ljj = (fin, puj / p) - ( pu ) (puj ) / p 

M tj = -^ ( kA 2 ) |S(v)| S'- (v) + (pA 2 |S(a)| S- (w)) A 
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( 21 ) 


( 22 ) 


(23) 


(24) 


The equation (20) yields six equations for one unknown coefficient Cd- The model coefficient is determined using a 
least square method 


( M .A) 


(25) 


In the present simulations, we used a simple top-hat filter with the filter-width two-times the grid size. This is 
equivalent of using k=V5. 


A. Solution Algorithm 

The governing equations were solved using higher order compact schemes. 28 29 Flow simulations were performed 
using an 8 th order compact scheme in the streamwise and spanwise directions and a 6 th order compact scheme in the 
wall normal direction. The orders are reduced to 4 th and 3 rd orders near the walls. We used higher order implicit 
filtering 30 of 10 th and 8 th orders in the homogeneous directions and in the normal directions, respectively to remove 
the high wavenumber contents in the solution. We employed a 3 rd -order total-variation-diminishing (TVD) Runge- 
Kutta scheme for time integration. We used a body-fitted curvilinear grid system in all of the simulations. The 
equations are transformed from the physical coordinate system (x, y, z) to the computational curvilinear coordinate 
system (£, q, Q. The grids were uniform in the streamwise and spanwise directions and were stretched in the normal 
directions close to the walls. The flow is maintained by applying a body fore e fi(t) in the streamwise direction. This 
is a function of time only and is determined at every step of the time marching by requiring that the average mass- 
flux remains constant, Eq. (26). We adapted the procedure described in Ziefle et al. 18 to determine the body force 
fi(t) at every time step. 


d 

dt 


J pudv = 0 

vol 


(26) 


III. Results 

A. DNS 

We first present the DNS results at a low Reynolds number of Rei,= 2800. We used (513, 257, 289) points in the 
streamwise, wall normal, and spanwise directions, respectively. The grid spacing in viscous units change along the 
streamwise direction due to the change in friction velocity along the wall. Figure 2(a) shows the variations of the 
grid spacing in wall units in the x, y , and z directions along the bottom wall and the grid spacing in the y direction in 
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the midway points between the bottom and top walls of the channel. The maximum grid spacing in the x- and z- 
directions are below 3 in the region up to x ~ 7.5 and increase to Ax + = 7.5 and Az + = 6.7 above the downstream hill 
due to a large increase in the friction velocity in this region. The minimum grid spacing at the wall varies between 
Ay + = 0.1 to 0.5 and the maximum spacing in the middle of the channel is Ay + = 3.3. Figure 2(b) shows the variation 
of the ratio of grid spacing in each direction to the Kolmogorov length scale, r\, along the normal direction at x = 
1.0. The maximum grid size is about 2.7 times the Kolmogorov scale. Hence these grid distributions are sufficient to 
resolve the flow field up to the dissipation scales. 

(a) (b) 




Figure 2. (a) Grid spacing in wall units along the bottom wall, and (b) ratio of grid spacing to Kolmogorov scale along the 

normal direction at x = 1.0. Re b = 2800. 


Figure 3 displays the instantaneous w-velocity field in the (x-y) plane at the spanwise location z = 0. The blue 
color denotes the reverse flow. The flow field can be divided into two regions. The first is the region below the line 
connecting the two crests of the hill and the second is the region above this line and below the top wall. The flow 
field in the second region appears as a plug flow through a two-dimensional channel. The first region is dominated 
by the reverse flow extending from one crest to the other. It is also seen that patches of attached flow exist inside 
and between the reverse flow regions. 



Figure 3. Instantaneous M-velocity Held in the (x-y) plane at z = 0. 


Figure 4(a) displays the instantaneous M-velocity field in the cross sectional (y-z ) plane across the separation 
bubble at the streamwise location x = 2.0. It is seen that the distorted shear layer is confined between the upper plug 
flow and the lower recirculation flow. Figure 4(b) shows the instantaneous M-velocity in the plan view along a fixed 
normal grid (j=5, y + ~ 2 ) plane ((x-z) plane) which is located a short distance above the lower wall. The white 
region displays the region with the negative velocity. The figure shows that except near the windward side of the 
downstream hill the streamwise velocity is negative in most of the domain. This phenomenon is common in 
separated flows 14 . The instantaneous separation regions extend to larger domains compared to those derived from 
the time averaged flow field. Figure 5(a) shows the variation of the instantaneous streamwise, u, and spanwise, w, 
velocities in the x-direction along a fixed normal grid line (J= 5) at a station z = 0. Similarly, Fig. 5(b) shows the 
velocities in the z-direction at the station x = 2.0 and y = 1.0. This station is located in the middle of the shear layer. 
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Figure 4. zi-velocity contours (a) in the cross sectional (y-z) plane at x = 2.0, and (b) in the plan view (x-z) plane at a fixed 

normal grid j= 5. 


(a) 



Figure 5. Variations of the instantaneous u, w velocity fields (a) in the x direction along a fixed grid line (j=5) at the station 

2 = 0 , and (b) in the z direction atx = 2.0 andy=1.0. 


The statistical quantities are obtained by averaging the solution in the spanwise direction and in time. The 
averaging in time was performed for about 40 flows through periods. Figure 6 shows the mean streamwise velocity 
contours and the streamlines from the current DNS. We observe backflow and massive separation behind the first 
hill, as expected. Figures 7(a) and 7(b) show the skin friction c/ and the pressure coefficient c p along the lower and 
the upper surfaces. We also include the DNS results 17 ' 19 obtained for the incompressible flow. In the current 
compressible computations at M 0 = 0.2, the flow separates at x sep = 0.23 and reattaches at x reatt = 5.45. The separation 
point location is close to the DNS 19 and the LES simulations 18 predicted value of 0.21. The current predicted 
reattachment point location is close to that predicted by the DNS 19 and the LES 18 simulations, where the 
reattachment occurred at x rmlt =5.4 and 5.3, respectively. The agreement of the skin friction coefficient and the 
pressure coefficient from the present simulations with the DNS 19 are quite good. The skin friction coefficient along 
the upper surface is higher compared to the incompressible DNS 19 results. The lower surface skin friction coefficient 
becomes negative after the separation and reaches a minimum value near the middle of the separation bubble. 
Beyond this minimum point, it slowly increases and takes a dip at the foot of the hill. The skin friction increases to 
larger values over the windward part of the hill and decreases at the crest. The pressure first increases and remains 
flat up to x ~ 2, and increases strongly up to the foot of the windward part of the hill due to the strong deceleration of 
the flow, x ~ 7.5, before it decreases to the inflow value. 
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Figure 6. Hill flow contours of the mean streamwise velocity and the streamlines, Re b =2800. 

(a) Skin friction (b) Pressure 




Figure 7. Hill flow variation of mean (a) skin friction coefficient, and (b) pressure coefficient along the lower and upper 

walls. Re b = 2800. 


In Figs. 8-10, the computed mean velocity <U> and the Favre-averaged Reynolds stresses <pu "u "> and <pu "v"> 
profiles are compared with the results of Breuer et al. 19 at five stations x = 0, 1, 2, 4, and 8. Mean velocity profiles 
and Reynolds stress profiles show excellent agreement with the DNS results of Breuer et al. 19 . This comparison 
validates the accuracy of the present computations. Figures 1 l(a-f) show all the Reynolds stress components profiles 
at six stations x = 0, 1, 2, 4, 5 and 8. We also plotted the mean velocity profiles in these figures. The first station x = 
0 is located at the hill crest and the flow separates slightly downstream of this point. The stations x = 1, 2, 4 and 5 
are situated inside the separation zone and the last station x = 8 is located on the windward side of the downstream 
hill. The first observation is that the Reynolds stresses in the separated region are confined to the separated shear 
layer region and they become comparatively small near the wall. The maximum values for the normal Reynolds 
stresses in the x-, y-, and z-directions at the station x = 1 are <pii "u "> =0.09, <pv "v "> =0.04, <pw "w "> =0.05 and the 
maximum Reynolds shear stress is <pu"v"> =0.03. These maximum values occur near the middle of the shear layer 
close to y ~ 1.0. The ratios of the three normal stresses are 1: 1/2.25 : 1/1.8. This shows that the normal stresses in the 
spanwise and the normal directions have almost the same magnitude. These ratios in a plane channel flow 31 are 1: 
1/9 : 1/6.25. This implies that the turbulence in separated shear layers is more isotropic than in wall bounded flows. It 
is also seen that when the flow approaches the reattachment point, the normal stresses in the streamwise and 
spanwise directions become almost the same near the wall. This suggests that the flow near the reattachment region 
behaves similarly to a stagnation flow impinging against a wall at an angle. The profiles near the downstream hill 
reveal that the maximum fluctuations occur in the spanwise direction compared to x- and y-directions. This was also 
observed in the simulations of Frohlich et al. 15 at a high Reynolds number of Re/, = 10595. This is due to the 
impinging of downstream moving eddies against the windward side of the downstream hill. This also generates 
substantial Reynolds stresses in the outer region above the crest as evident in Figs. 1 1 (a) and (f). 
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(a) x=0 


(b) x=1 


(c) x=2 


(d) x=4 


(e) x=8 




Figure 8. Hill flow mean streamwise velocity profiles at different stations. Re b = 2800. 

(a) x=0 (b) x=1 (c) x=2 (d) x=4 (e) x=8 
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Figure 9. Hill flow streamwise normal Reynolds stress profiles at different stations. Re b = 2800. 



(a) x=0 


(b) x— 1 


(c) x=2 


(d) x=4 


(e) x=8 
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(a) x=0 






The transport equations for the Reynolds stresses and the turbulent kinetic energy are given in Ref. 32. The 
equation for the turbulent kinetic energy can be written as 


1 


K = —OU; U, 

2 # ii 

d JL,Y-D- e -C 

dt 


(27) 


The left hand side of Eq. (27) is the local rate of change of turbulent kinetic energy and the terms on the right 
hand side are: (1) production, (2) diffusion, (3) dissipation, and (4) advection. In the statistical steady state 
conditions, the local rate of change should be zero and these four terms should balance each other to yield zero sum. 
A non-zero left hand side points to lack of grid resolutions in the simulations. This non-zero value is termed as 
numerical dissipation and in turbulent simulations this balance is added to the physical dissipation and considered as 
the total physical dissipation in the analysis. Figures 12(a-e) depict the magnitudes of the four terms on the right 
hand side of this equation and the balance or the total of these four terms at five stations x = 0, 1, 2, 4 and 8. It is 
seen that the balance is below 3% of the dissipation in the whole domain. This confirms that the grid used in this 
DNS resolved all the scales including the dissipation scales at this low Reynolds number. The first observation is 
that turbulent production is the maximum at all the stations. The production peaks near the separation point x = 0 
and 1 to a value of 0.10 and decreases to 0.05, 0.03 and 0.02 at the stations x = 2, 4 and 8, respectively. It is also 
observed that these peaks occur in the middle of the shear layer. Let us consider the balance of the four terms at the 
station x = 1, which is downstream of the separation point x = 0.23. The production peaks at a value of 0.10 neary ~ 
1. This is balanced in the negative side by the advection, diffusion and dissipation terms. The magnitudes of these 
terms are 0.01, 0.06 and 0.03, respectively. Hence diffusion is about two times larger than the dissipation term. 
When we move towards the lower wall, all the terms first decrease to zero and at the wall dissipation again 
increases. The dissipation at the wall is balanced by the positive diffusion term. This description remains the same 
up to the separation point except that the diffusion term decreases faster than the dissipation term. 
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(a) x = 0 


(b) x = 1 


(c) x = 2 





(d) x = 4 


(e) x = 8 




Figure 12. Balance of different terms in the turbulent kinetic equation. Reb = 2800, DNS. 

We also plotted the contours of the average statistical quantities: (a) production, (b) dissipation, (c) ratio of 
production to dissipation, (d) kinetic energy, and (e) shear stress in Figs. 13(a-e). These quantities play important 
roles in RANS modeling and development. The production is maximum near the start of the separated shear layer as 
observed earlier and is confined to a region along the separation line. The maximum production occurs near x = 0.6 
and y = 1.06, which is downstream of the separation point and slightly above the hill crest. Beyond the reattachment 
point, the production is concentrated in the middle of the channel at a height of y ~ 1. There is negative production 
near the reattachment region and along the shoulder of the hill. It is also noted that there exists a high production 
region above the downstream hill. The maximum dissipation occurs along the wall on the windward side of the hill. 
Outside the wall region, dissipation is confined to the shear layer and to its proximity. The ratio of production to 
dissipation is shown in Fig. 13(c). This ratio takes a peak value of 4.5 near the start of the separated shear layer at x 
= 0.6 and y = 1.06. It decreases to a value around 2 near x - 2 and v ~ 1. Recall that this ratio is about 1.8 in the 
buffer region of a turbulent channel flow. It is also observed that this ratio remains at a constant value of about 1.5 
up to the second hill along the line y ~ 1. Figure 13(d) shows the turbulent kinetic energy distribution in the 
separated turbulent boundary layer. The kinetic energy is concentrated along the shear layer and it peaks near.v ~1.4 
and v - 1.06. After the reattachment point the energy is confined to region near y - 1. This agrees with the earlier 
observations from the turbulent production contours. After the reattachment point, turbulence is confined to a region 
away from and parallel to the wall at a height ofy - 1. The turbulent shear stress contours in Fig. 13(e) also convey 
the same picture that turbulent shear stress is concentrated near the start of the shear layer and near the y ~ 1 region 
atx- 1.4. 
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Figure 13. Hill flow contours of the turbulent statistical quantities (a) Production, (b) Dissipation, (c) 
Production/Dissipation, (d) Kinetic energy, and (e) Shear stress. Re b = 2800. 
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B. LES 

We have implemented the dynamic Smagorinsky model (DSM) into the existing code and performed the LES 
simulations for the two Reynolds numbers Re b = 2800 and 10595 using the higher order compact scheme. The grids 
that were used in the DNS and the LES are given in the Table 1. 

Table 1. Grid sizes used in the DNS and LES simulations 


Re 

DNS 

LES 

2800 

513*257*289 

Grid 1 129*129*73 

Grid 2 257*257*145 

10595 


Grid 1 193*193*145 

Grid 2 257*257*217 


(1) Case 1: Re b = 2800 

Figure 14 shows the variations of the grid spacing in wall units for the two grids, Grid 1 and Grid 2, in the x, y, 
and z directions along the bottom wall and the grid spacing in the y direction in the midway points between the 
bottom and top walls of the channel. The maximum grid spacing in the x- and z-directions for Grid 1 are below 12 in 
the region up to x ~ 7.5 and increase to Ax + = 30 and Az + = 25 above the downstream hill due to large increases in the 
friction velocity in this region. The minimum grid spacing at the wall varies between Ay + = 0.70 to 2.0 and the 
maximum spacing in the middle of the channel is Ay + = 6. Similarly, the maximum spacing in the x- and z-directions 
for Grid 2 are below 6 in the region up to x ~ 7.5. The minimum grid spacing at the wall varies between Ay + = 0.1 to 
0.9. These values are finer than the recommended values 32 of Ay + < 2, Ax + ~ 50-150 and zlz + ~15-40 for wall 
resolved LES. Especially, the values for Grid 2 are much finer than the recommended values. 



Figure 14. Grid spacing in wall units along the bottom wall. 

Figures 15(a-d) show the contours and the streamlines obtained from (a) DNS, (b) LES without the SGS model 
(ILES(l)) using Grid 1, (c) LES with the dynamic Smagorinsky model (LES(l)) using Grid 1, and (d) LES with the 
dynamic Smagorinsky model (LES(2)) using Grid 2. We did not include the results obtained without the SGS model 
using Grid 2 because the results are almost the same as LES(2). Qualitatively, the large scale features such as 
separation point, reattachment point, bubble size and shape, and the velocity contours look almost the same in all 
four figures. The bubble looks slightly smaller and the reattachment occurs earlier in the results obtained with 
ILES(l). Figure 16 depicts variation of the skin friction at the lower wall. The reattachment point is located at x reatt = 
5.4, 5.0, 5.2 and 5.4 from the four simulations DNS, ILES(l), LES(l) and LES(2), respectively. The corresponding 
separation points are located close to x ~ 0.23, 0.23, 0.24 and 0.25, respectively. The implicit LES simulation using 
the coarser grid (ILES(l)) predicts an early reattachment point compared to other cases. This is due to the fact that at 
this grid level the dissipation scales are not resolved and the energy dissipated by the numerical scheme is not 
sufficient to keep the turbulent kinetic energy at the correct level. Flence, there exist large turbulent fluctuations in 
the shear layer and this strong mixing causes the shear layer to reattach earlier. In these simulations, we employed 
uniform grid distributions in the x- and z-directions. In earlier simulations, 17 ’ 18 ’ 19 a stretched grid was used in the x- 
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direction. This makes the grid finer near the separation point and coarser in the middle part of the channel. Whether 
the stretched grid distribution causes changes in separation and reattachment points will be evaluated in the future. 
Figure 17 shows the variation of the sub-grid eddy viscosity at different streamwise stations x = 0.0, 1.0, 2.0, 4.0, 
8.0. As expected, the eddy viscosity is large in the separated shear layer region. However, the computed maximum 
eddy viscosity is only 50% of the molecular viscosity for the Grid 1 and is about 30% for the Grid 2. This may be 
the reason that at this low Reynolds number the SGS model did not affect the dynamics of the turbulence. 


Y 


(a) DNS 



(b) LES No Model (Grid 1) 



(c) LES with Dynamic Model (Grid 1) 



(d) LES with Dynamic Model (Grid 2) 



Figure 15. Hill flow contours of the mean streamwise velocity and the streamlines, Re b =2800. (a) DNS, (b) ILES (Grid 1), 

(c) LES with DSM (Grid 1), and (d) LES with DSM (Grid 2). 
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Figure 16. Hill flow variation of mean skin friction coefficient along the lower wall. Re b = 2800. 



Figure 17. Variation of the sub-grid eddy viscosity at different stations. 

Figures (18-21) depict the comparison of the mean and turbulent quantities obtained from the four simulations at 
different streamwise locations x=0, 1, 2, 4. The mean streamwise and normal velocities obtained from the ILES and 
LES simulations are in excellent agreement with the DNS results at x = 0, 1 and 2. The velocities at x = 4 are slightly 
smaller than the DNS results near the wall region. The Reynolds stress components also show good agreement with 
the DNS results. The maximum normal Reynolds stress <pu "u "> at x = 4 obtained with LES(l) is under predicted 
by about 5% compared to the DNS results. 
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(a) 


(*>) 




Figure 18. Comparison of mean streamwise velocity <U> profiles obtained from the DNS, ILES and LES simulations at 

the stations (a) x = 0 and 1, and (b) x = 2 and 4. 


(a) (b) 




Figure 19. Comparison of mean vertical velocity <V> profiles obtained from the DNS, ILES and LES simulations at the 

stations (a) x = 0 and 1, and (b) x = 2 and 4. 
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(a) 




Figure 20. Comparison of turbulent intensity < pii "u "> profiles obtained from the DNS, ILES and LES simulations at the 

stations (a) x = 0, and (b) x = 4. 


(a) (b) 




Figure 21. Comparison of turbulent shear stress < pii "v"> profiles obtained from the DNS, ILES and LES simulations at 

the stations (a) x = 0, and (b) x = 4. 


(2) Case 2: Ret, = 10595 

Figures (22-33) similarly show the LES results at the higher Reynolds number of 10595. Figure 22 shows the 
variations of the grid spacing in wall units for the two grids, Grid 1 and Grid 2. The maximum grid spacing in the x- 
and e-directions for Grid 1 are below 20 in the region up to x ~ 7.5 and increase to Ax + = 55 and Az + = 36 above the 
downstream hill. The minimum grid spacing at the wall varies between Ay + = 0.20 to 4.0 and the maximum spacing 
in the middle of the channel is Ay + = 10. Similarly, the maximum spacing in the x- and z-directions for Grid 2 are 
below 15 in the region up to x ~ 7.5. The minimum grid spacing at the wall varies between Ay + = 0.20 to 1.0. These 
values for Grid 2 are finer than the recommended values 31 in the whole domain. The grid spacing near the wall for 
Grid 1 is slightly higher than the recommended values. 
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50 



Figure 22. Grid spacing in wall units along the bottom wall, Re = 10595. 

Figures 23(a-d) show the contours and the streamlines obtained from (a) LES without the SGS model (ILES(l)) 
using Grid 1, (b) LES(l) with the dynamic Smagorinsky model using Grid 1, (c) LES without the SGS model 
(ILES(2)) using Grid 2 and (d) LES(2) with the dynamic Smagorinsky model using Grid 2. The separation bubble 
appears slightly smaller and the reattachment point occurs earlier in the results obtained with ILES(l) compared to 
LES(l). Figure 24(a) depicts the variation of the skin friction at the lower wall. Figure 24(b) shows the pressure 
coefficient along the lower and upper walls. The pressure coefficients agree very well with the Breuer et al. 19 data 
except the implicit LES simulations with Grid 1. The pressure coefficient up to x ~ 2 is lower than that was 
predicted by other simulations. The reattachmcnt point is located at x reatt = 3.75, 4.2, 4.4 and 4.4 from the four 
simulations ILES(l), LES(l), ILES(2) and LES(2), respectively. The separation points are located close to x ~ 0.22, 
0.23, 0.22 and 0.20 from the four simulations. The separation point location is close to the DNS 19 and the LES 17 ' 18 
predicted value of 0.21. The predicted reattachmcnt point locations by the DNS 17 ' 19 and the LES 18 vary between 4.6- 
4.7. The implicit LES simulation using the coarser grid predicts an early reattachmcnt point compared to other cases. 
This was also observed in the simulations at the lower Reynolds number Reb = 2800. As we discussed earlier, this is 
due to the low molecular dissipation of the turbulent kinetic energy by the resolved scales. Figure 25 shows the 
variation of the sub-grid eddy viscosity at different streamwise stations x = 0.0, 1.0, 2.0, 4.0, and 8.0. As we 
observed earlier, the eddy viscosity is large in the separated shear layer region. The computed maximum eddy 
viscosity is about 1.8 times the molecular viscosity for the coarser Grid 1 and it is about 1.2 for Grid 2. Flence, with 
the increasing Reynolds number, the effect of the SGS increases and this is reflected in the separation lengths 
obtained with and without the SGS model. 
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3.0 


(a) LES No Model Grid 1 



(c) LES No Model Grid 2 




Figure 23. Hill flow contours of the mean streamwise velocity and the streamlines, Re=10595. (a) ILES(l), (b) LES(l), (c) 

ILES(2) and (c) LES(2). 
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(a) Skin friction (b) Pressure 




Figure 24. Hill flow variation of mean (a) skin friction coefficient, and (b) pressure coefficient along the lower and upper 

walls. Re = 10595. 



Figure 25. Variation of the sub-grid eddy viscosity at different stations. Re = 10595. 

Figures (26-29) depict the comparison of the mean and turbulent quantities obtained from the four simulations, 
ILES(l), LES(l), ILES(2), and LES(2) at different streamwise locations x = 0 and 4 with the data of Breuer et al. 19 
The mean velocity profiles show good agreement with the LES results of Ref. 19. The results from the dynamic LES 
are better than that obtained with the ILES(l). The normal mean velocity at x = 4 is slightly under predicted 
compared to the data from Ref. 19. The Reynolds stresses obtained with the dynamic LES also depict good 
agreement at these stations with the data from Ref. 19. 
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(a) X = 0 


(b) X = 4 




Figure 26. Comparison of mean streamwise velocity <U> profiles obtained from the ILES(l), LES(l), ILES(2) and LES(2) 

simulations at the stations (a) x = 0, and (b) x = 4. 


(a) X = 0 


(b) X = 4 




Figure 27. Comparison of mean vertical velocity <V> profiles obtained from the ILES(l), LES(l), ILES(2) and LES(2) 

simulations at the stations (a) x = 0 , and (b) x = 4. 
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(a) X = 0 


(b) X = 4 




Figure 28. Comparison of turbulent intensity <pu"u"> profiles obtained from the ILES(l), LES(l), ILES(2) and LES(2) 

simulations at the stations (a) x = 0, and (b) x = 4. 


(a) X = 0 <b) X = 4 




Figure 29. Comparison of turbulent shear stress <pu"v"> profiles obtained from the ILES(l), LES(l), ILES(2) and 

LES(2) simulations at the stations (a) x = 0, and (b) x = 4. 

Figures 30(a-f) show all the Reynolds stress components profiles at six stations x = 0, 1, 2, 4, 5 and 8 obtained 
from LES(2). We also plotted the mean velocity profiles in these figures. The first station x = 0 is located at the hill 
crest and the flow separates slightly downstream of this point. The stations x = 1, 2, and 4 are situated inside the 
separation zone, the station x = 5 is located downstream of the reattachment point, and the last station x = 8 is 
located on the windward side of the downstream hill. The conclusions are similar to what we observed in the low 
Reynolds number case (Figs, ll(a-f)). The maximum Reynolds stresses in the separated region are confined to the 
separated shear layer region and they become small near the wall. The maximum values for the normal Reynolds 
stresses in the x-, y-, and z-directions at the station x = 1 are < pu "u "> = 0.08, < pv "v "> = 0.05, < pw "w "> = 0.06 
and the maximum Reynolds shear stress is < pu"v"> = 0.04. These maximum values occur near the middle of the 
shear layer close to y ~ 1.0. The ratios of the three normal stresses are 1: 1/1.6 : 1/1 .3. As we observed earlier, the 
turbulence in shear layers is more isotropic than in wall bounded flows. It is also noted by comparing these ratios 
with that for the low Reynolds number Ret,= 2800 case that the flow in the shear layer region becomes more 
isotropic with the increasing Reynolds number. It is also seen as we observed for the low Reynolds number case that 


22 


when the flow approaches the reattachment point. Figs. 30(d) and (e), the normal stresses in the streamwise and 
spanwise directions become almost the same near the wall. As we described in the low Reynolds number case, this 
is due to the impinging of the shear layer against the bottom wall. The profiles near the downstream hill reveal that 
the maximum fluctuations occur in the spanwise direction compared to x- and indirections. These observations agree 
with the simulations of Frohlich et al. 15 



Figure 30. Normal and shear Reynolds stress components profiles at different stations. Re b = 10595. 

Figures 31(a-e) depict the magnitudes of the four terms in the energy equation Eq. (27), and the balance or the 
total of these four terms at five stations x = 0, 1, 2, 4 and 8 obtained from LES(2). The observations are again similar 
to what we discussed in the low Reynolds number case, Figs. 12(a-e). However, in this high Reynolds number case, 
it is observed that the balance term is not zero for the grid we used. This points to the lack of grid resolution in the 
simulations. The maximum deficit occurs in the separated shear layer regions at x = 1 and 2. The percentages of the 
maximum deficit compared to the resolved dissipation are 25, 60, 50, 45, 40 at the stations x = 0, 1, 2, 4, and 8, 
respectively. These percentage differences were also observed in other simulations. 14 ' I5 ' 19 The maximum production 
at the hill crest location x = 0 is about 0.2, which is two times higher than that in the low Reynolds number case. 
However, downstream of the separation point, the magnitudes and the profiles for the four terms are almost the same 
as that in the low Reynolds number case. The maximum production is about 0.10 at the station x = 1 and it decreases 
to 0.04, 0.02, and 0.01 at the stations x = 2, 4, and 8, respectively. It is also observed that these peaks occur in the 
middle of the shear layer. Comparing the figures for the low and high Reynolds numbers, Fig. 12 and 31, it is seen 
that the role of different terms in the turbulent kinetic equation budget is almost the same in both cases. Near the 
separation point region, x= 0 and 1, the production is balanced by the advection, diffusion, and dissipation terms. It 
is also seen that diffusion is about two times larger than the dissipation term in this region. When we move towards 
the lower wall, all the terms first decrease to zero and at the wall dissipation again increases and this is balanced by 
the positive diffusion term. 

We also again plotted the contours of the average statistical quantities: (a) turbulent production, (b) turbulent 
dissipation, (c) ratio of turbulent production to dissipation, (d) turbulent kinetic energy, and (e) turbulent shear stress 
in Figs. 32(a-e). The observations are similar to the low Reynolds number case. The production is maximum near 
the start of the separated shear layer and is confined to a region along the separation line. Beyond the reattachment 
point x ~ 4.4, the production is concentrated in the middle of the channel at a height of y ~ 1. There is negative 
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production near the reattachment region and along the shoulder of the hill. The maximum dissipation occurs along 
the wall on the windward side of the hill. Outside the wall region, dissipation is confined to the shear layer and to its 
proximity. The ratio of production to dissipation is shown in Fig. 32(c). This ratio takes a peak value of 6.5 near the 
start of the separated shear layer x ~ 0.4 and y ~ 1.0. It decreases to a value around 2 near .r ~ 6 and y ~ 1. Figure 
32(d) shows the turbulent kinetic energy distribution in the separated turbulent boundary layer. The kinetic energy is 
concentrated along the shear layer and it peaks near x ~ 1 and y ~ 1.0. After the reattachment point, the energy is 
confined to a region near y ~1. After the reattachment point, turbulence is confined to a region away from and 
parallel to the wall at a height of v ~ 1. The turbulent shear stress contours in Fig. 32(e) also convey the same picture 
that turbulent shear stress is concentrated near the start of the shear layer and near the y ~ 1 region at the foot of the 
hill. 

To investigate the effects of the SGS in the LES simulations, we compared the Reynolds stresses and the 
different terms in the turbulent kinetic energy equation obtained with LES with no model ILES(l) and LES with the 
dynamic Smagorinsky model LES(l) using Grid 1. Figures 33(a) and (b) show the comparison at the station x = 1. 
The solid and dashed lines denote the results for ILES(l) and LES(l), respectively. We also plotted the production 
term obtained from LES(2) as the thick black line in Fig. 33(b). It is seen that turbulent production obtained with 
ILES(l) is 30% higher compared to LES(l) and LES(2) at this station. The second observation is that the resolved 
dissipation with ILES(l) is two times smaller than that obtained with LES(l). This will manifest into higher 
Reynolds stresses in the ILES(l) computations compared to LES(l) computations. Figure 33(a) shows that all the 
Reynolds stresses are higher in the ILES(l) compared to LES(l). This confirms the role and the importance of the 
SGS in the coarse grid LES simulations. The sub-grid stresses extract the energy from the large scales and dissipate 
it. 


(a) x = 0 (b) x = 1 (c) x = 2 





(d) x = 4 (e) x = 8 




Figure 31. Balance of different terms in the turbulent kinetic equation. Re b = 10595. 
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Figure 32. Hill flow contours of the turbulent statistical quantities (a) Production, (b) Dissipation, (c) 
Production/Dissipation, (d) Kinetic energy, and (e) Shear stress. Re b = 10595. 
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(a) X=1 Reynolds stresses (b) X = 1 K.E Budget 



Figure 33. Comparison of (a) Reynolds stresses and (b) balance of different terms in the turbulent kinetic equation 
obtained from ILES(l) (solid lines), LES(l) (dashed lines). Black thick solid line in Fig. (b) production only from LES(2). 

IV. Conclusions 

We investigated the turbulent flow through a channel flow with a constriction (2-D hill) using DNS and LES 
calculations at two Reynolds numbers Ret = 2800 and 10595. The Navier-Stokes equations are solved using a higher 
order compact scheme for space discretization and a 3 rd order Runge-Kutta scheme for time discretization. We 
employed the dynamic Smagorinsky sub-grid scale model in the LES calculations. DNS simulations were performed 
using a fine grid at the low Reynolds number of 2800 and the results are compared with the reference data of Breuer 
et al. 19 The separation and reattachment points, skin friction, pressure coefficient, mean flow, and the turbulent 
quantities agree well with the reference data. The instantaneous flow fields identified two distinct regions. One is the 
region between the line joining the crests and the bottom wall and the other is above this line and the top wall. The 
upper region resembles a plug flow through a channel. The bottom region is dominated by the recirculation flow and 
large eddies. These two regions are separated by a strong shear layer. The results showed that the Reynolds stresses 
peak in the separated shear layer region near the separation point. The simulations also showed that the normal 
Reynolds stresses in the wall normal direction and in the spanwise directions have almost the same magnitude in the 
shear layer region. Close to the reattachment point the magnitudes of the normal stresses in the streamwise and the 
spanwise directions become almost the same. The budget of the terms in the turbulent kinetic energy equation 
showed that the production is balanced by the diffusion, dissipation and advection terms in the shear layer region 
near the separation point. It is also observed that the diffusion is two times larger than the dissipation in the 
proximity of the separation point. The ratio of production to dissipation takes a maximum value of 4.5 near the 
separation point. Another general observation is that all the quantities, shear stress, turbulent kinetic energy, 
production, remain at a reasonable level in a region away from and parallel to the wall at the hill height. 

LES simulations were performed at two Reynolds numbers Re* = 2800 and 10595 using two grids in each case. 
We also performed simulations with and without the SGS model using the same grid. The mean flow and the 
turbulent quantities agree well with the reference data. The predicted reattachment point is about 4% smaller 
compared to the reference data 19 . The computed eddy viscosity showed that for the grids and the algorithm used in 
these simulations, the effect of SGS is small at the lowest Reynolds number of 2800. The effect of the SGS increases 
with increasing Reynolds number. The Reynolds stresses and the budget terms have similar behavior at both low 
and high Reynolds numbers. The ratio of production to dissipation increases to 6.0 in the high Reynolds number 
case. The simulations with and without the SGS model showed that the production is larger without the model 
compared to with the model. The dissipation decreases without the model compared to with model. These two 
effects cause the Reynolds stresses to increase and force the flow to reattach earlier compared to that with the model. 
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